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Abstract. 

The UUersma model for the damped harmonic oscillator is coupled to the 
quantised electromagnetic field. All material parameters and interaction strengths 
are allowed to depend on position. The ensuing Hamiltonian is expressed 
in terms of canonical fields, and diagonalised by performing a normal-mode 
expansion. The commutation relations of the diagonalising operators are in 
agreement with the canonical commutation relations. For the proof we replace 
all sums of normal modes by complex integrals with the help of the residue 
theorem. The same technique helps us to explicitly calculate the quantum 
evolution of all canonical and electromagnetic fields. We identify the dielectric 
constant and the Green function of the wave equation for the electric field. Both 
functions are meromorphic in the complex frequency plane. The solution of the 
extended UUersma model is in keeping with well-known phenomenological rules for 
setting up quantum electrodynamics in an absorptive and spatially inhomogeneous 
dielectric. To establish this fundamental justification, we subject the reservoir of 
independent harmonic oscillators to a continuum limit. The resonant frequencies 
of the reservoir are smeared out over the real axis. Consequently, the poles of 
both the dielectric constant and the Green function unite to form a branch cut. 
Performing an analytic continuation beyond this branch cut, we find that the 
long-time behaviour of the quantised electric field is completely determined by the 
sources of the reservoir. Through a Riemann-Lebesgue argument we demonstrate 
that the field itself tends to zero, whereas its quantum fluctuations stay alive. We 
argue that the last feature may have important consequences for application of 
entanglement and related processes in quantum devices. 
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1. Introduction 

More than a decade ago, Huttner and Barnett ^ published a valuable contribution 
to the complicated and long-established subject of quantum electrodynamics in 
nonrclativistic macroscopic matter. Opting for a canonical setting, they carried 
out a fundamental derivation of the quantised electromagnetic field in an absorptive 
dielectric medium. To open up the possibility of working analytically, they avoided 
making direct contact with the atomic level. Instead, they described the properties of 
the dielectric with the help of harmonic oscillators. The influence of absorption was 
mimicked by coupling the dielectric to a continuum of harmonic oscillators. These 
reasonable simplifications gave rise to a so-called damped-polariton model that could 
be exactly solved. The quadratic Hamiltonian was diagonalised by means of Fano's 
method |H] and Fourier transformation. 

The solution for the quantised electromagnetic field in an absorptive dielectric 
was welcomed by a large community. Phenomenological quantisation schemes, the 
practical value of which was beyond dispute, could be put on a solid microscopic 
foundation. This possibility marked the beginning of a period of progress, in which 
our knowledge of QED in or near macroscopic matter was considerably broadened. For 
a lot of different optical media and experimental geometries, ranging from magnetic 
materials to beam splitters, quantisation of the electromagnetic field was successfully 
carried out. At present, the macroscopic formulation of QED is well established. A 
few years ago, it was comprehensively reviewed 

Over the years, it was pointed out several times |7] that the justification of 
Huttner and Barnett has one major shortcoming. The dielectric is assumed to be 
homogeneous in space, whereas in many experimental situations the electromagnetic 
field is substantially influenced by spatial inhomogeneities. One thus would like to 
extend the solution of the homogeneous damped-polariton model to the case in which 
all model parameters depend on position. This is a far from trivial enterprise, because 
the lack of translational invariance deprives us of the possibility to perform Fourier 
transformation. On the other hand, the quadratic character of the Hamiltonian 
remains intact, so there are no mathematical indications for a failure of Fano's 
procedure when passing over to the inhomogeneous case. Indeed, in two companion 
papers it is demonstrated that the inhomogeneous damped-polariton model can be 
solved as well. One can employ either Fano's procedure or an alternative method 
that is based on Laplace transformation 9 . 

The extension of the work of Huttner and Barnett to dielectrics with spatial 
inhomogeneities does not complete the program of underpinning phenomenological 
quantisation rules. This judgement, which is the main motivation for undertaking 
the present study, is supported by several arguments. First, standard works 
on conservative QED and classical fields IT^ insist on reducing 

the collection of canonical variables to an enumerable lot before commencing a 
canonical quantisation procedure. In doing so for the inhomogeneous case, we are 
automatically led to a complete set of normal-mode functions, on the basis of which 
the electromagnetic fields can be expanded. We thus identify the functions that replace 
the plane waves of the homogeneous case. Moreover, we find the natural decomposition 
of the Green function IHUHI belonging to the wave equation for the electric field. 

Our second reservation with regard to the damped-polariton model concerns 
the loss mechanism. Instead of employing a continuum of oscillators right from 
the beginning, we should postpone the transition to the continuum as long as 
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possible. This allows for a scrutiny of the mathematical origin of the absorptive 
behaviour. Also, we can completely clarify the relation between irreversibility of the 
dynamics and causality of the dielectric constant. Our third argument pertains to a 
technical observation on Fano's method in the presence of continua. In solving for the 
diagonalising operators one must introduce a formal contribution that is proportional 
to a Dirac delta function. We wish to improve upon this approach. To that end, 
we must refrain from employing distributions Pj. Fourth and last, we remark that 
already back in the 1960s quantum dissipation was studied with the help of harmonic 
oscillators nSl~CH|- It is important to find out whether the older oscillator models 
corroborate the predictions of the damped-polariton model. 

We can meet the four suggestions advanced above by exchanging the continuum of 
the damped-polariton model for a reservoir containing a finite collection of independent 
harmonic oscillators. In essence, we add an electromagnetic sector to the Ullersma 
model |T7] for a damped harmonic oscillator. As expected, we can still diagonalise 
the Hamiltonian of the extended Ullersma model. Eventually, we come up with a 
mathematical limit that restores the continuum. 

In section 2 we specify the Hamiltonian of our model. In preparation of a swift 
canonical quantisation procedure, we solve the classical evolution equations for the 
canonical fields by invoking a normal-mode expansion. Quantisation of the dynamics 
happens in section 3. With each normal mode we associate a quantised harmonic 
oscillator. Next, making use of complex integration, we determine the evolution in 
time of all canonical fields. This permits us to identify the dielectric constant as 
well as the Green function. Subsequently, we turn on absorption by subjecting the 
reservoir to a continuum limit. We witness how poles unite and radically change the 
analytic structure of both the dielectric constant and the Green function. We derive 
the inhomogeneous counterpart of the solution that was obtained by Huttner and 
Barnett. Much attention is devoted to the long-time behaviour of the electric field. 
Section 4 contains a summary and guides the reader to the main results of our paper. 

We close this introduction with some technical remarks. In performing partial 
integrations for canonical fields and related quantities, we tacitly assume that there 
are no contributions from the boundaries of the dielectric. Our aim to keep the 
treatment free from distributions forces us to utilise a considerable amount of basic 
function theory; all of the accompanying calculus is transferred to two appendices. In 
an attempt to present clear formulas, we omit spatial arguments whenever possible. 
The subscript L{T) denotes that one should take the longitudinal (transverse) part of 
a vector field or tensor field. Last, we make use of rationalised mks units throughout. 

2. Classical treatment 

In this section we derive the classical Hamiltonian of our model. We let us be guided 
by the Lagrange formalism, so our treatment bears a canonical character from the very 
outset. As usual, Hamilton's equations furnish the evolution laws for the canonical 
fields. In solving these, we rely on a normal-mode expansion. It appears that the 
dynamics is governed by an eigenvalue problem for a self-adjoint differential operator 
in three dimensions. For a spatially homogeneous dielectric the eigenvectors can be 
constructed from plane waves. To verify that the normal modes neatly decouple, we 
compute the Hamiltonian on the basis of the solutions for the canonical fields. 
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2.1. The model 

The Ullersma model ^7] describes the interaction between a single harmonic oscillator 
and a reservoir made up by an array of N independent harmonic oscillators. To fulfil 
our purposes, we assign to all N+1 oscillators a spatial dependence. The corresponding 
position vector r covers a finite volume V . As a further extension of the Ullersma 
configuration we introduce an electromagnetic sector, denoting the electric field and 
magnetic field at time t as E(r, i) and B(r,t), respectively. In vacuum these fields 
give rise to the standard electromagnetic Lagrangean density 

CEM^\e,,^^~^B\ (I) 

The presence of a dielectric medium of volume V is taken into account by the privileged 
harmonic oscillator of the extended Ullersma model. It has mass density p(r) and 
frequency CL'o(r). The field Qo(r, t) measures its displacement. The Lagrangean density 
of the dielectric reads 

- IpQI ~ Ipu^oQl (2) 

After suitable scaling, all harmonic oscillators of the reservoir have a mass density of 
p(r) as well. Their frequencies and displacement fields are equal to cc'„(r) and Q„(r, t), 
respectively. From ((2J) we see that 

^ /I 1 \ 

= E 2^^" " o'^'^nQl (3) 

ri=l ^ ^ 

is the Lagrangean density of the free reservoir. 

In the electric-dipole approximation the electric field induces in the dielectric 
a polarisation density V = — aQo, where a is positive. In a strictly microscopic 
theory — a(r) would be a local electronic charge density. The interaction between 
the electromagnetic fields and the dielectric yields a contribution of —as4> + js • A to 
the Lagrangean density. In absence of free charges and currents, one may substitute 
CTs = — V-P and js = V for the sources. The electromagnetic potentials are determined 
by the definitions E = — V(/) — A and B = V x A, supplemented with the choice 
Ai = 0, which is equivalent to the Coulomb gauge. 

As long as it allows for energy exchange, the precise form of the interaction 
between the dielectric and the reservoir is not of physical interest to us. Introducing 
a coupling /3n(r), we link the displacement field Qo to the time derivative Q„ rather 
than the field Q„ itself. This departure from the Ullersma model is common in 
studies of dissipative QED 1 . The two interactions of our model make the following 
contribution to the Lagrangean density: 

N 

Ci = a(V0) • Qo - aA • Qo - 51 /^"Qo • Qn- (4) 

n=l 

We have added the total derivative V • (a^Qo ) , so that the nabla operator acts on the 
scalar potential. 

As the Lagrangean L — J dr{CEM + + + Ci) does not depend on 0, 
one immediately finds \J(j) — — (ckQo)L/£o from the corresponding Euler-Lagrange 
equation. Owing to the Coulomb gauge, the character of the first term of Q is 
purely electrostatic. The remaining Euler-Lagrange equations provide us with the 
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inhomogeneous Maxwell equation for the vector potential, and the equations of motion 
for all A'' + 1 displacement fields. 
Defining canonical momenta as 

5L „ 5L 



n 



5A (5Q„ 
for < n < iV, we obtain the Hamiltonian 

N 

H ' ' ' 



j dr |^A-n + ^^Q„-P„^ 



as 



H = dr 



N 



2eo 



2^0 



2p 



2eo 



2p 



(5) 



(6) 



— + (V X A)2 + — + — (aQo)i + -Pg + -pQ^Ql 



(7) 



For the sum ujq + J2n=i Pnl abbreviation Wg will be in use. The Hamiltonian 
of our model being available, we can start investigating the evolution of the canonical 
fields. 



2.2. Solution oj Hamilton's equations 

Since we work in a finite volume, we may try to unravel the dynamics with the help 
of an enumerable set of independent normal modes. We shall see that each mode has 
a specific spatial structure and oscillates at a specific frequency. For that reason, the 
modes must be labelled by both a spatial index k and an integer which enumerates 
the mode frequencies i7(fc, I) for k fixed. The index j identifies the canonical fields 
with the auxiliary fields Zj in the following manner: 

{n. A, Po, Qo, P„, Q„} = {Zi, Z2, Z3, Z4, Z } (8) 

with n = 1, 2, 3, . . . , A^. Now we can put forward our normal-mode expansion 

(r, t)^Y^ c{k, l)aj {k,l; r)e-*"('='')* + cc. (9) 

k,l 

The frequencies n{k,l) are positive by definition. The coefficients c{k,l) could 
be absorbed in the mode amplitudes a.j, but this is inconvenient in view of the 
quantisation procedure lying ahead of us. Because of the Coulomb gauge the 
amplitudes ai and 3.2 are transverse vector fields. 

Upon substituting @ into Hamilton's equations and carefully evaluating 
functional derivatives, one arrives at 

iSiai = Aa2 + — aa + -ag 

Mo \ P J T \P 

i^lsL2 = ai 

eo 



N 



lila^i = pwoa4 H (q!£14)^+_> — a5„ 

Co — \ P 
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is2a4 — a2 a3 

P P 
iila^n = P^^n^en 

i^aen = - — a4 a^n (10) 

P P 

with 1 < 71 < A^. The arguments of a^ and are identical to those appearing in Q. 
We emphasize that the model parameters a, /?„, p, wq, and cl;„ depend on position. 

We set out to derive a wave equation for the amplitude e{k,l;r) of the electric 
field. The expansion 

E(r, t)=Yl Oe(fc, I; Oe-^"^'^^')* + cc (11) 

should match with the definition of E in terms of the electromagnetic potentials. This 
brings us to the prescription 

eoe = -ai + (Q;a4)^ . (f2) 

The above relation enables us to present the solution of the algebraic part of (|10|l as 

i 

ai = -eoBT a2 = -t^gt 



ia iafl 
as = TrGT - TTT^e a4 



a 



The new function h{s) implicitly depends on position, and is given by 

In the next section, the real-valued argument s will be replaced by a complex variable. 
A swift verification of p3l) may take place through substitution into (|10|l . Indeed all 
five algebraic equations are satisfied if (|12|l is used. 

The first equation of pU|) . containing the Laplacian operator, has not been 
considered as yet. After substitution of (|13|l and employment of (|12|l it reduces to 

c^V X (V X e) = n'^e{n)e. (15) 

This is the standard wave equation for the electric field, which appears in 
electromagnetic theory. The function 

e{s) = 1 - -4n (16) 

eoph{s) 

plays the role of dielectric constant. Its spatial dependence is made explicit in the 
following. 

In H15(l a linear differential operator is at work. It is given by 

L{s; r)v(r) = -c^V x [V x v(r)] + s'^e{s; r)v(r). (17) 

The parameter s is arbitrary, but still real-valued. The vector field v belongs to a 
Hilbert space of complex and square integrable functions, with scalar product defined 
as 

(vi,V2) - J dr vJ(r)-V2(r). (18) 



Oscillator model for dissipative QED in an inhomogeneous dielectric 



7 



The integration is confined to the finite volume V . Since L{s) is self-adjoint, its 
eigenvalues \{k, s) are real. Second, its eigenvectors u(A:, s) possess the orthonormality 
property 

{u{k,s)Mk\s))=5kk'. (19) 

Obviously, H19|) is not affected by degeneracy. Eigenspaces of higher dimension should 
be subjected to an orthonormalisation procedure. Third, the eigenvectors make up a 
basis for the Hilbert space. We have a decomposition of the tensorial delta function 
at our disposal, given by 

Y,n*{k,s-v')n{k,s-v)^d{v' -v). (20) 

k 

The reality of the left-hand side follows from the invariance of the tensorial delta 
function under interchange of its arguments and indices. Last, as L{s) is quadratic in 
the parameter s, the eigenvectors for s and —s can be chosen identical to each other, 
so that the relation 

u(/c, — s; r) = u(fc, s; r) (21) 

may be assumed. 

Of course, for s equal to a mode frequency f2, the eigenvalue problem 

L(s; r)u(fc, s; r) = A(fc, s)u(fc, s; r) (22) 

should be in keeping with the wave equation (|15|) . We therefore make the choice 

e(fc, /; r) = w{k, l)u{k, n{k, I); r). (23) 

The weight w will be evaluated in due course. We furthermore require that the 
eigenvalue A(fc, s) be equal to zero if s coincides with a mode frequency. In other 
words, by solving the equation 

A(fc,s)=0 (24) 

with s real-valued, we gather all mode frequencies. From (|19|l and (|22|) one infers 

A(fc, s) = (u(fc, s),s^e{s)u{k, s)) -c^{\/ x u(fc, s), V x u{k, s)). (25) 

A partial integration has been performed. In appendix A we use (|25ll to take a closer 
look at the solutions of (|24|l . From (|21|) we see that A(fc, s) is even in s, so only pairs 
s = ±ri occur. The negative solutions generate the cc terms of 

In summary, Hamilton's equations have been solved by means of a normal-mode 
expansion. The spatial structure of the normal modes is described by the eigenvalue 
problem (|22|l . whereas the mode frequencies are specified by the constraint (|24|l . Now 
we should keep our promise to demonstrate that the modes do not interact with each 
other. 



2. 3. Computation of the Hamiltonian 

We insert the normal-mode expansion (O for all canonical fields into the Hamiltonian 
(0. Subsequently, we employ H13|) and H23|) to convert all amplitudes a.j{k,l) into 
eigenvectors u{k, ft). The arguments (fc, I) of are omitted. Likewise, by il' the mode 
frequency Q{k',l') is meant. If we perform a partial integration, the wave equation 
H15|l allows us to get rid of all differential operators. The longitudinal contributions 
of JTI) can be eliminated via the relation 

[e(O)u(fc,17)]i = (26) 
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which is a consequence of ((T3|) and if^ . 

The Hamiltonian is equal to a sum over indices fc, Z, k' , I'. We can interchange 
these to make all summands maximally symmetric. Then two classes of summands 
remain: those oscillating rapidly at frequency 17 + f2', and those oscillating slowly at 
frequency fl ~ fl' . For each summand of the first class there is a complex conjugate, 
which may be ignored in the sequel. For clarity we remark that summands oscillating 
at frequency — f2 + f2' can be transferred to the second class. A simple interchange of 
summation indices suffices. 

The summands of rapid oscillation will be treated first. All of these can be 
expressed in terms of h functions, owing to the relation 

P'iuj'^^Lul - SS') sh{s) + s'h{s') 2/2 / -2 

„^lP'('^n-s')('^n-s'2) S + S' 

The sum on the left-hand side originates from Q . The dummies s and s' are set equal 
to and ri', respectively. Verification of (|27|l is straightforward after substitution of 
the definition (|14|l on the right-hand side. If we eliminate all h functions in favour of 
dielectric functions , the summands of rapid oscillation yield a form that vanishes 
after use of the identity 

(u*(fc', n'), n'^e{n')u{k, n) ~ n^e{n)uik, n)) = o. (28) 

One proves (|28|l by taking the scalar product of u*(fc', ft') and the wave equation H15() . 
Once again, partial integration is indispensable. 

We have found that only the summands of slow oscillation contribute to the 
Hamiltonian. As long as differs from Q', these summands can be handled in the 
same manner as discussed above. There is only one difference. In order to finalise the 
calculation, one needs the identity 

(u(fc', n'), n'^e{n')u{k, n) - n'^e{n)u{k, n)) = o (29) 

instead of The proofs of and follow the same path. Altogether, for 

J7 ^ ri' the summands of slow oscillation add up to zero as well. 
The case f2 = 17' calls for employment of the relation 

+ _dshis) (30) 



Now elimination of h in favour of e no longer produces a null result. We are led to 

„2 



^1 P^(W^-S2)2 ds 

of h in favour of e no longe 
H^2eo J2 'cik,l)c*{k',l')wik,l)w*ik',l')(uik',s) /^y^^''K {k,s)\ . (31) 

Because of the constraint 17 = O' the summation carries a prime. As the wave equation 
l|15|l is only valid for discrete values of 17, it may not be differentiated with respect to 
this variable. Therefore, further manipulations of (|31|l must take place in a meticulous 
manner. From partial integration and H15|) we learn 



, ds^e{s) , 
u{k ,s), — p— u(fc,s^ 



^{u{k',s),L{s)u{k,s)) 



(32) 



s=Sl 



with the condition 17 = 17' in force. The eigenvalue equation (|22|l holds true for all real 
values of s, so it may be utilised on the right-hand side of (|32l) . Bearing in mind that 
the eigenvectors are orthonormal, and moreover, that at fixed k the correspondence 
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between r2(fc, I) and I is one-to-one, we can finish our computation of the Hamihonian. 
The diagonal form 

d\{k, s) 



i7-2eo5]|c(fc,0l'k(fc,0l' 



k,l 



(33) 



s=n{k,i) 



is found. In accordance with our expectations, the normal modes do not interact with 
each other. 

We choose the weights as 

-1/2 



w{k, I) 



-1/2 



n{k,i) 



dX{k, s) 
ds 



s=n{k,i} 



so that 133() becomes 

H = J2^{k,lMk,l)\' 



(34) 



(35) 



k,l 



This is the Hamiltonian of an enumerable collection of independent harmonic 
oscillators. The quantity c{k, l) is sometimes called the amplitude of the normal mode 
k,l ^5]. With the normal-mode expansion © and the Hamiltonian (|35|l in hand, we 
fully understand how the classical dynamics of our model works. 



2.4- Homogeneous dielectric 

Analytic solution of the eigenvalue problem (|22|l will be impossible, except for special 
cases. One of these is a spatially homogeneous dielectric. The model parameters 
no longer depend on position. One can dispose of the parameter s figuring in the 
eigenvectors of (I22f) . If one does not care for rotational or other spatial symmetries, 
then a simple solution for the eigenvectors is provided by plane waves. 

A cube of side L serves as the volume V of the dielectric. The spatial index k 
splits up into a polarisation index = 1,2,3 and a wave vector q = 27rm/L, with 
rrij any integer. The eigenvectors u(fc;r) are given by L^'^/^o^ cxp(zq • r), where 63 
equals q and the set {61, 62, 63} is orthonormal. In order to fulfil (|20|l the continuum 
limit L ^ 00 should be taken. 

The solution of l|22|) and l|25|) is composed of a longitudinal part and a transverse 
part. The longitudinal and transverse mode frequencies obey the equations 

e{n) ^ Q^s{Q) = c2q2. (36) 

The fact that the dielectric is invariant under rotations causes the degeneracy in H36f) . 



3. Quantum treatment 



In the previous section we recognised that the classical evolution of our extended 
UUersma model is controlled by an enumerable set of independent harmonic oscillators. 
For the quantisation of each oscillator we resort to the old method of Dirac. The 
advantage is that the quantised Hamiltonian is delivered to us in diagonal shape. 
Since we work in a canonical setting throughout, we must ascertain that the canonical 
operators obey the canonical commutation relations. Completion of this job teaches 
us that discrete sums over normal modes can be transformed into complex integrals, if 
in the eigenvalue problem 1)22(1 the parameter s is replaced by a complex variable. The 
transition to complex integration is the key to disclosing the quantum evolution of the 
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extended UUersma model. For all Heisenberg operators, including the electric field 
and the electric displacement, integral representations are obtained. The dynamics is 
in agreement with Heisenberg's equations. To bring about dissipation in the dielectric, 
we make use of a continuum limit, which allots to the reservoir an uncountable number 
of degrees of freedom. Once the energy sink is activated, we can identify the causal 
dielectric function. 

3.1. Dirac quantisation 

From here onwards, all capitals refer to quantummechanical operators, the Green 
function G, the operator L, and the frequency O excepted. Following Dirac, we 
associate with each normal mode k,l of energy ^l{k,l)\c{k,l)\'^ a harmonic oscillator. 
We define ladder operators through 

c(fc, I) h^'^C{k, I) c*{k, I) h^'^C\k, I) (37) 

and quantise by postulating 

[C(k, I), C\k\ I')] = Skk'Sw [C{k, I), C{k\ I')] = 0. (38) 

Upon symmetrizing properly, we obtain the quantised counterpart of the Hamiltonian 
as 

H = ^"^VL{k,l) [C\k,l)C{k,l) + C{k,l)C\k,l)\ . (39) 

Each eigenvalue Vl is positive. From and (|37|) we deduce that all canonical operators 
are represented by the expansion 

(r, t) = n^/^ "^&j{k,l- r)C(fc, Oe-'^e^'')* + he. (40) 

k,l 

The amplitudes a.j are completely determined by the results established in the previous 
section. In calculating the quantised Hamiltonian, one may also depart from Q. After 
substitution of 140() into I?)), one enters a rather lengthy road. The manipulations are 
essentially the same as in section 2.3, so there is no need for any comments. One 
indeed retrieves the diagonal form H39(l . 

3.2. Canonical commutation relations 

On our way to the swift quantisation procedure 138(1 . we let us be guided by the 
Lagrange formalism. Therefore, one may rightfully ask whether all is well with the 
canonical commutation relations. These read 

[A(r',t),n(r,t)] =iMT(r'-r) 

[Qo(r',i),Po(r,i)] - *M(r' - r) 

[Q„Kr',i),P„(r,i)] = zM„'„<5(r' - r). (41) 

All other commutators of canonical fields equal zero. In the following we shall only 
verify the upper two commutators (|41|l . For the other 19 cases new problems do not 
arise. Incidentally, for some commutators verification is trivial. 

In virtue of ((38|l and ((40(1 the upper two conditions ((41(1 can be brought onto the 

form 

E T^^''T^^^ ^' "^^Tik, I; r) + cc = drir' - r) (42) 
k,i ^^^"'^f 



Oscillator model for dissipative QED in an inhomogeneous dielectric 



11 



Ea[r)a[r')e*{k, I; r') 
p{v')n{k,i)Hnik,iy,r') 

The relations and suggest a recourse to the residue theorem. In moving to 
the complex plane we observe the rule 

u(fc, s; r) u(fc, z; r) u*(fc, s; r) — > [u(fc, z*; r)]* (44) 

where u(fc, z; r) is the solution of H21|) and 122|) with the replacement s — > z carried 
out. In appendix A we argue that in (|44|l two entire functions of z figure. Moreover, 
we make plausible that the function A^^(fc,z) is meromorphic, and that its poles are 
given by z = ±n{k, I), for k fixed. 

Upon performing in H42|l and l|43|) the transition to complex integration, we find 
the following three sufficient conditions: 

E/ -T7f^[uT(A,z*;r')]*UT(fc,z;r)=^T(r'-r) (45) 
J, Jci ^I'^J ^) 

The contour Ci is composed of a set of circles running in counterclockwise sense (see 
figure Bl in appendix B). The Ith circle encloses the pole on the positive real axis 
at z = fl{k,l), with k fixed. Everywhere else in the interior of Ci the integrands 
are analytic, because the circles can be chosen as small as we like. Hence, use of the 
residue theorem in (|151l - llT7|l indeed reproduces the sums over mode frequencies that 
are contained in H42I) and 143|l . 

The conditions H45|l -147 () can be proved by performing a series of contour 
deformations. We defer this purely technical exercise to appendix B. Right now the 
reader should stay focused on the transition to complex integration, as practised above. 
We plan to exploit that skill in computing the time evolution of the canonical fields 
and other operators of physical interest. This is the goal of the next subsection. 

3. 3. Solution of the extended Ullersma model 

To get a full picture of the dynamics of our quantum system, we have to specify the 
evolution of any set of initial canonical operators. In short, we have to establish the 
mapping between the times t = Q and t for all canonical operators. To that end, 
we observe that any quantum operator can be written as a linear combination of the 
canonical operators at time zero. We apply this statement to C''{k,l), and use (|40|l 
as well as (|41|l to make all coeflticients explicit. We are led to the expansion 

ai(fc, I; r) • A(r, 0) - a2(fc, I; r) • n(r, 0) + a3(fc, I; r) • Qo(r, 0) 



n^{k,l)e{k,l;r) 
h{n{k,l)-r) 



exik^l; 



cc = ^(r' - r). (43) 



N N 

-a4(fc, I; r) • Po(r, 0) + ^ a5„(fc, I; r) • Q„(r, 0) - ^ a6„(fc, l; r) • P„(r, 0) 



(48) 
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One may wonder whether the representation (|48|l comphes with the quantisation 
prescription H38|l . By invoking (|13|) . (|23|) . and H41|) . one derives a set of consistency 
relations that is equivalent to (|28|l . (|32|l . and (|34|l . 

After substitution of (|48|) into ()40|l and use of the normahsation H34|l . we can 
express the vector potential as 

The source vector must be constructed from the initial canonical fields. One has 
c2j(s; r) = ~*6osA(r, 0) + n(r, 0) + [a(r)Qo(r, 0)]^ - ^'"WQ"('-'0) 

\(50) 



isQ;(r)Po(r, 0) ^ sQ;(r)/?„(r) 
p(r)/i(s;r) p{r)h{s;Y)[s^ - uHr)] 



zo;2(r)Q„(r,0)-^P„(r,0) 
p(r) 



The intermediate result H49|l paves the way for the residue theorem, in a similar vein 
as before. The following result is reached: 

A(r, t)= j dr' -^e— *Gt(z; r, r') • j(z; r'). (51) 

The operation T refers to the argument r. The contour C3 encloses the real axis by 
means of two straight lines running from +00 + irj to —00 + ir], and from —00 — irj to 
+00 — ir], where r] is infinitcsimally positive. From H22I) we see that the Green function, 
defined as 

G(z; r, r') - ^ ^(F^"^^' ^' ""^ I"^^' ""'^l* ^^^^ 
satisfies the partial differential equation 

c-2l(z; r)G(z; r, r') = ^(r-r') (53) 

where z must lie on the contour C3. As shown in appendix B, the integrand of (|51|l 
has the same analytic structure as X~^{k,z). From (|14|l it is clear that in H5U|) the 
factor of [z^ — a;^(r)]~^ does not give rise to any poles. Hence, in (|51|l the contour Ci 
could be exchanged for C3 without paying a price. 

For the displacement field Qo the discrete solution is given by 

OJrt)- V / , , f c^g«(r)e— *u(fc,5;r)u*(fc,g;rO-j(g;rO ) +v,r(^A) 
Qo(r, t)- l^Jdr^ ^^^^^ ^^^^^ + he. (54) 

The transition to complex integration is immediate, provided that Ci is chosen as 
contour. In deforming the latter to C3, attention must be paid to denominators in 
which a double factor of h figures. Correction terms show up, which can be found by 
eliminating the nasty denominators with the help of the auxiliary function f , defined 
in appendix B. One should proceed along the same lines as for the proof of H47|) . This 
results in 

Qo(r,<) = - / dr' [ -^f^W£2!!G(z;r,r')-j(z;r') 

-/ :^^T{c'j(2;r) + zeo2A(r,0)-n(r,0)- [a(r)Qo(r,0)]^}. (55) 
2tti za[r) 
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The local character of the correction term stems from apphcation of the completeness 
relation (j^OJl. 

The displacement field of the reservoir contains two local correction terms. The 
solution comes out as 

^ / N /",//" dz z^a(r)Bn(r)e^^^* ^, r. ., 
Q„ r,i = - / dr' / J- ^^"V z;r,r' -j z;r' 

- / V- (^ 7\l 2 27^{a(^;r)+^eo^A(r,0)-^(r,0)~Kr)Qo(r,0)]J,} 

Jc^2TTa{r)p{r)[z^ -uji{r)] ' 

dz 



c,2^z [z^-Lol{v)] 



zQ„(r,0) + ^P„(r,0) 



(56) 



where one has n — 1,2,3, ...,A^. The last correction term covers the special case of 
— Pn — 0. The fields A, Qo, and Q„ provide us with the following solutions for the 
canonical momenta: 

n = eoA Po = pQo - aA P„ = pQ„ - /3„Qo (57) 

with n = 1, 2, 3, . . . , A'^. These expressions are obtained from the definitions |SJ- 

All solutions for the canonical operators reproduce the initial condition if the 
choice t = is made. This can be verified on the basis of the identities underlying 
the canonical commutation relations, such as (|45|l and H47|l . The transformation 
z — > —z* shows that the condition of self-adjointness is respected as well. One needs 
the symmetry property 

[G(-z*;r,r')]* = G(z;r,r') (58) 

where z belongs to C3. For the proof one combines (|52(l with the reciprocity relation 

G(z;r',r),, = G(z;r,r')»j. (59) 

This last result follows by taking the scalar product of H53|) with G(z;r, r"), and 
performing a partial integration. 

All canonical fields obey the Heisenberg equation ifiiij = [Z^, H]. To demonstrate 
this, one first calculates the commutator with H on the basis of the canonical 
commutation relations. Subsequently, one employs (|57|l to eliminate all canonical 
momenta. Last, one substitutes the integral solutions 151|) . (|55|l . and (|56|l . Contour 
deformations are not required; use of the partial differential equation H53|l is sufficient. 

To compute the time evolution of the electric field we repeat the program 
described above for Hll() . This brings us to 

E(r, t) = - I dv' I -^ze-''%{z- r, r') • j(z; r'). (60) 
Combination of (|55|l and (|60|) gives for the electric displacement 

D(r,t) = eoE(r,t)-a(r)Qo(r,<)=D(r,G)+ / — c2j(z;r) 

-fdv'f ^z£(z;r)e-^^*G(z;r,r')-j(z;r'). (61) 

From H6()|l and (|61(l we conclude that one may indeed regard e(z; r) as a dielectric 
function. With the help of (|17|l one can clarify the status of the local contribution 
to the electric displacement. It guarantees that the field DL(r,f) does not exist, as 
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prescribed by Maxwell's equations. Indeed, if one takes the divergence of (|51|l and 
inserts (|53|l . the right-hand side reduces to the divergence of D(r, 0), which equals 
zero. 

In the electric field and electric displacement were calculated for the case of 
a spatially homogeneous dielectric, coupled to an uncountable number of harmonic 
oscillators. We can extend these results to the case of an inhomogeneous dielectric 
by considering and H61|l for a reservoir, the eigenfrequencies of which make up a 
dense set. This will be done in the next subsection. 

5.^. Continuum limit 

For finite N the solution of the extended UUersma model describes reversible dynamics. 
Therefore, the energy exchange between dielectric and reservoir is everlasting. On the 
other hand, for many experiments on optical properties of dielectrics, absorption of 
photons is omnipresent. Hence, there is still a gap between the solutions of the previous 
subsection and experiment. Our dielectric medium does exhibit damping phenomena 
if we manage to create an irreversible energy flow into the reservoir. For that purpose, 
the reservoir should possess an uncountable number of degrees of freedom. Such a 
continuum comes into existence upon defining 



for n = 1, 2, 3, . . . , iV, and taking the limit A, — > oo. Since n becomes arbitrarily 
large, the ratio n/A may be treated as a continuous variable w. Note that the spatial 
dependence of (jJ„ does not leave any traces. The subscript c indicates that the 
continuum limit is taken. 

We assume that for all nonnegative to the function /3(w) is regular and smooth. 
If it decays faster than 1/v^ for large lo, then the definition below {T)) implies that 
the continuum limit of cjq exists. The analytic properties of the function 



radically differ from those of h{z), given in (|14|l . The zeros of h{z) have united so as 
to generate a branch cut on the real axis. The same mechanism is witnessed for the 
dielectric function ec{z) and the eigenvalue \c{k^z). In appendix A we demonstrate 
that all of the afore-mentioned functions are analytic and nonzero outside the real 
axis, that is to say, both on the branch in the upper half plane and on the branch in 
the lower half plane. We assume that the analytic properties of the eigenvector n{k, z) 
are not affected by the continuum limit. 

Now we are well prepared to find out how the solution of the extended UUersma 
model, which was derived in the previous subsection, behaves under the continuum 
limit. We focus on the electric field, specified in (|60|l . Treatment of the electric 
displacement and other fields goes by the same methodology. As mentioned earlier, the 
contour C3 is composed of two straight lines enclosing the real axis. We parametrise C3 
as / dzf{z) — dLu[f{Lu — irj) — fito + irj)] , where 77 is infinitesimally positive. Next, 
we simplify the integrand with the help of the symmetry [G(z*; r, r')]* — G{z;r,r'). 
Application of rule (|62|l then leads to 



Ai/2Q„(r,t) = Q(n/A;r,i) 
Ai/2/3„(r)=/?(n/A;r) 



Ai/2p„(r,t)=P(n/A;r,i) 
ujn{r) =n/A 



(62) 




(63) 
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dv' 



duj 



*Ci;a(r')Im 



+ 



X Im 



duj 
Treo 

Gc{i^ + ii]\ r, r') 



Gc(w + i77;r, r') 



hc{uj + iri]v') 



■in 



f , du; r , 

dr / duj 

J J-oo ^'^0 Jq 



p{r') 



hc^LU + irj; r')[(LU + ir])^ — lu 



12] 



■ iR{uj,uj';r'). 



(64) 



The continuum Green function 
replacements z ^ lu + irj and A 

by 

c^j£;Af('^;r) = 



Gc{i^ + iri;r,r') is determined by H52|) . with the 
-> Ac carried out. The new source vectors are given 

- zeowA(r, 0) + n(r, 0) + [a(r)Qo(r, 0)]^ 



c^jD(w;r) =iuj' 



'Qo(r,0) 



UJ 



■Po(r,0) 



c^iniuj, Lo'; r) = iu'^CI{lo'- r, 0) - ^P(^'; r, 0). 



(65) 



In the above solution is derived on the basis of Laplace transformation. In that 
paper, the reservoir contains a continuum of oscillators right from the beginning. 

By taking the continuum limit as demonstrated above, we fully extend the results 
of to the case of an inhomogeneous dielectric. However, there exists an alternative 
manner to implement the continuum limit. Instead of keeping the upper and lower 
part of C3 together, one can decide to sever these parts from each other. The ensuing 
representations are useful to analyse how fields behave for long times. We keep the 
time strictly positive and focus again on the electric field. For some integrands figuring 
in it is wise to perform first the substitution z/[f{z)\ = [z^ — f {z)] / [z f {z)] + l/z, 
where f{z) stands for Ac(fc, z) or hc{z). One then isolates the term l/z, which decays 
slowly as |z| becomes large. The convergence of the remaining integrand improves. 

The lower part of C3 can be closed by means of the arc z = i?exp(i0), with 
— TT < (p < 0. Owing to the choice t > 0, the exponential factor exp{—izt) causes 
the integrand to shrink to zero on the arc. In the interior of the closed contour the 
integrand is analytic, so the lower part of C3 does not make any contribution. The 
upper part of C3 yields an integral over all real frequencies uj, which is convergent in 
virtue of the above substitution. Of course, for the term l/z one cannot treat the 
upper and lower parts of C3 separately; instead, the residue theorem offers a way out. 

If we execute the above instructions for the electric field, then the following 
expression emerges: 



E(r,t) -E(r,0) 



dr' 



-»c^ta(i")[^c(t^ + i?7;r) 



duj 



{lo + iri)hc{LU + irj; r) 



■Qo(r,0) 



dr' 



-00 27rzeo 

w2a(r')Qo(r',0) 



wGc(w + irj\ r, r') 



IT] 



-5{r^r') 



C^hc{LO 

duj _,„, 



i77;r') 



2a(r') 



2'Kie, 







■^Po(r',0) 



p{r')hc{ijJ + ir]] 
r') 



-Gc(w + ir7;r,r') 



did' 



P{lo' 



[UJ ■ 



iriY 



—iR{uj,uj';r') . (66) 



Oscillator model for dissipative QED in an inhomogeneous dielectric 



16 



As before, rj is infinitesimally positive. In computing hc{Lo + irj] r) the upper branch of 
the function H63|l must be used. Consequently, each integrand of H66|) is analytic in the 
upper half of the complex w plane. For t = one may then replace the contour by an 
arc u) — Rexp{i(j>), with < (j) < tt. On account of (|14|l . (|20|l and H22|l . each integrand 
decays as or faster, so there are no integrals surviving the choice t — 0. In the 

continuum limit the initial condition for the electric field is still satisfied. 

To uncover irreversible behaviour, we have to analyse H66|l for large times. If the 
symmetry relation 

[P^{-z*)]* =P'{z) (67) 

holds true, and P{z) is well-behaved at the origin, then the analytic continuation of 
ec{uj + iri) and Ac(fc, w + zry) can be effectuated. Thus the infinitesimal increment ir] of 
the argument uj becomes redundant. From p6(l and the continuation of hc{uj + irj) we 
obtain the symmetry relation [£c(— w*)]* = Sc{uj), which is in accordance with classical 
theory of the dielectric constant jEl- All this is shown in appendix A. There we also 
prove that edto) is analytic and different from zero for all Imuj > 0. Incidentally, if 
the symmetry (|67|) is absent, indeed circumstances exist under which the process of 
analytic continuation fails l2Uj . 

Now everything is ready to shift in (|66|) the contour into the lower half plane. The 
poles at Lo — —irj yield residues that erase the initial condition for the electric field. 
We adapt 7 such that all integrands accompanying A, 11, Qo, or Po, are analytic in 
the strip —7 < Imcj < 0. Then one may integrate along the line Imw = —7 instead of 
the real axis. The parameter 7 surely differs from zero, otherwise the above process 
of analytic continuation would fail. The exponential exp(— iwi) produces a factor of 
exp(— 7i), which induces absorption in the dielectric, as desired. By comparing (|6()|l 
and H66|l one can trace the origin of the irreversible behaviour. It resides in the fact 
that, as a result of the continuum limit, the analytic properties of h{z) undergo a 
radical change. 

Use of (|16|l in H66|l gives a dielectric function ed'jj) that is analytic in the upper 
half plane. Therefore, irreversibility for t — > 00 is coupled to causality for the dielectric 
function. Indeed, if one takes the continuum limit for t negative, one has to operate 
on the lower branch of (|63|l . Then the dielectric function is analytic in the lower 
half of the complex uj plane. One encounters the combination of irreversibility for 
t — !■ —00 and anti-causality for the dielectric function. Of course, these conclusions 
do not depend on 1)67(1 . In absence of this symmetry one may resort to arguments of 
Riemann-Lebesgue type. 

The terms of (|66f) that contain the canonical fields of the reservoir have not yet 
been considered. The corresponding integrand has poles at a; = ±0;' — irj, the residues 
of which generate oscillating contributions. Upon using the symmetry relations for Gc 
and /ic, we arrive at the following asymptotic expression for the electric field, valid for 
large times: 

E(r,i)~ 

- / dv' dce-^-'^ I ""^^'^^"j;;^"' ''^^ G.(c.; r, ■ c.; + he (68) 

The phase (p — aigihc) results from eliminating the coupling parameters a and (3 in 
favour of the dielectric function. To do this, the usual distributional calculus should 
be applied to (|63|l . The imaginary part of the integrand on the right-hand side of (|63|l 
is proportional to a Dirac delta function. 
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In virtue of a Riemann-Lebesgue argument, the expectation value of the integral 
of Ht)8|) vanishes for large times. Hence, the same is true for the electric-field 
operator E(r, t). Still, the oscillatory contributions of the reservoir qualitatively differ 
from those exhibiting exponential damping. The oscillations bring about quantum 
fluctuations that do not die out for large times. To make this explicit, we model the 
initial correlations in the reservoir as 

(P(tj; r, 0)P(cj'; r', 0)) = \ {V^iuj: r, 0)) 5{uj - J)5{r - r') (69) 

where the last 8 function is a tensor. The brackets indicate that an expectation value 
is taken. For the auto-correlations of Q(w;r,0) the above model is assumed as well. 
Cross-correlations of P(ti-'; r, 0) and Q(w; r, 0) are discarded. From the electric field we 
construct a quadratic form, take the expectation value, and process the result with the 
help of and Utilising again a Riemann-Lebesgue argument and introducing 
the energy density of the reservoir as 

Hfl(a;;r) = ip-i(r)p2(a;; r, 0) -t- ia;V(r)Q'(w; r, 0) (70) 

we find the asymptotic result 
lini (E(r,i)E(r',t)) = 

t—^OQ 

dv" / -^yQ,{uo-ry) ■ G:(c.;r",r')Im[£,(a;;r")] {Ti r[u: ; v" ))+ cc. (71) 

In the continuum limit any expectation value of the electric field will eventually decay 
to zero, but the quantum fluctuations stay alive. They are fuelled by the energy that 
is available in the reservoir. 



4. Summary and conclusion 

The quantisation of the electromagnetic field in the presence of nonrelativistic 
macroscopic matter constitutes a problem of formidable magnitude. In principle, any 
serious treatment should start from the minimal-coupling Hamiltonian that takes into 
account all electromagnetic interactions at the atomic level. It needs no argument that 
such a rigorous approach is scarcely possible, unless one is prepared to deploy heavy 
numerical means as soon as it comes to predicting experimental findings. However, one 
then misses the opportunity to see what kind of mathematical mechanisms are at work 
behind such processes as spontaneous emission or scattering of photons. Therefore, 
already during the early years of the subject people were looking for shortcuts so 
as to obtain concise and transparent theories Nowadays a generally accepted 
quantisation scheme of phenomenological nature is available. 

Because of their experimental relevance, rules for QED in absorptive matter 
deserve to be put on a firm foundation. This idea was pursued by Huttner and Barnett 
^ . Sacrificing the direct contact with the atomic level, they solved a so-called damped- 
polariton model. Their microscopic expression for the quantised electromagnetic 
field in an absorptive dielectric initiated a lot of activity on the construction of 
phenomenological quantisation rules. Despite this progress, a for experimentalists 
most important issue is still pending, namely, the microscopic quantisation in a 
dielectric with both spatial inhomogeneities and losses. In addition to that, it is unclear 
how the damped-polariton model relates to the old oscillator models for dissipative 
quantum dynamics that were proposed in the 1960s |15|-|18|. 
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Motivated by the last remarks, we exchange in this work the continuum of the 
damped-polariton model for a reservoir that consists of a finite number of independent 
harmonic oscillators. We thus consider the UUersma configuration extended with 
an electromagnetic sector. For all parameters a spatial dependence is admitted, so 
that the oscillator model at hand is capable of describing spatial inhomogeneities. To 
give the reader the chance to make a direct comparison between the quantum and the 
classical evolution of the model, we first solve Hamilton's equations for the classical 
canonical fields. The latter are expressed as sums of independent normal modes. 
By associating with each normal mode a harmonic oscillator, we can invoke Dirac's 
method for the quantisation of the dynamics. The corresponding ladder operators can 
be employed to cast the Hamiltonian into diagonal form. Hence, the time evolution of 
all fields can be made explicit. Having derived the solution of the extended UUersma 
model, we can identify the inhomogeneous dielectric function. Also, we can meet one 
of our prime objectives, the generalisation of the results of Huttner and Barnett PP 
to the case of an inhomogeneous dielectric. This job is completed by subjecting our 
discrete reservoir to a continuum limit, which makes the coupling parameter of each 
separate oscillator infinitely weak, and at the same time converts the collection of 
eigenfrequencies into a dense set. 

The discrete character of the UUersma reservoir brings us important advantages. 
We can completely avoid the mathematical complexities that accompany the use of 
distributions. For example, we need not adapt the solution of the evolution equations 
through addition of an unknown term, containing a delta function El IE] • To derive 
the results of this paper only simple technical tools are required. We apply linear 
algebra to prove orthogonality and completeness for the eigenvectors of operator L, 
which governs the spatial structure of the normal modes. Hence, these eigenvectors are 
the natural candidates for composing the Green function, as becomes manifest when 
computing the solutions for the canonical fields. The other standard ingredient we call 
in, is the residue theorem. It allows us to reduce all discrete sums of normal modes 
to complex integrals. The ensuing contour Ci is rather awkward, as it is composed of 
separate circles around poles on the real axis. Fortunately, Ci gives way to frequency 
integrals upon performing contour deformation. In doing so, we generate correction 
terms of local nature, in which the Green function no longer figures. These terms can 
be interpreted as Langevin noise operators The local correction to the electric 
displacement ensures that the divergence of the latter equals zero, in accordance with 
Maxwell's equations. The easy access to the mathematical structure of our solutions 
enables us to scrutinise their behaviour under the continuum limit. The function h 
plays a central role. Its poles gradually cover the real axis, and finally give rise to 
a branch cut. The dielectric function is analytic on both sides of the cut, but the 
sign of the time t determines on which branch the fields differ from zero. This leads 
to the conclusion that (anti)causality of the dielectric function is intimately linked to 
irreversibility of the dynamics for t approaching (minus) infinity. 

A further advantage offered by the UUersma reservoir is the connection with 
recent work on decoherence and entanglement. Several papers describe how a reservoir 
of the UUersma type can bring about decoherence of quantum superpositions and 
entanglement between qubits 125-123 ■ Results on these processes are valuable with 
an eye to quantum computing and other future applications. Inevitably, in any real 
device electromagnetic fields will participate in the interactions. Therefore, it would 
be worthwhile to study the reservoir-induced quantum phenomena in the presence 
of an electromagnetic sector. One would like to assess the influence of the quantum 
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noise that is produced by the electromagnetic field. As we saw in this work, the energy 
density of the reservoir fuels the quantum fluctuations of the electric field. They outlive 
any decoherence or entanglement. The possible fragility of these quantum processes 
can surely be investigated with the aid of the extended UUersma model. Especially 
the influence of temperature must be critically regarded. 

Appendix A. Results for h and related functions 

In the main text we exploited some analytic properties of the functions h[z; r), e(z; r), 
and A(fc, z). For the first two functions, all properties can be proved by means of simple 
techniques, as we shall discuss below. When we turn to A(fc, z), however, the going 
becomes heavier, due to the fact that the partial differential equation (|22|) intervenes. 
We shall need a few plausibility arguments in order to make progress. A rigorous 
mathematical investigation of (|22|l lies outside the scope of this paper. 

The function h{z] r), defined in (|14|l . is equal to the ratio of two finite polynomials 
in z. Hence, h~^{z]Y) is a meromorphic function of z, the poles of which follow by 
solving the equation h{z; r) = 0. The imaginary part of h{z] r) is given by 



is manifest. We see that the square must be real, otherwise /i(z;r) cannot equal 
zero. For z^ <Q the inequality 



comes into play, where h is real. Clearly, h{z] r) can equal zero only if is real and 
positive. We confirm our surmise that h^^{z]v) is meromorphic, with all poles lying 
on the real axis, symmetric with respect to the origin. 

The function /ic(z;r), the continuum counterpart of /i(z;r), is specified in H63() . 
As long as Imz differs from zero, the factor of (w^ — z^) cannot render the denominator 
of (|53|) zero. This guarantees that hc{z\ r) is analytic outside the real axis, i.e., on each 
of its two branches. Since hc{z] r) is obtained by applying a limiting process to /i(z; r), 
the inequalities HA.2|I and (|A.3p are valid for /ic(z;r) as well. Therefore, /ic(z;r) does 
not vanish outside the real axis. 

In virtue of (|16|l the dielectric function ec{z\ r) inherits the branch cut of /ic(z; r), 
and is analytic on each of its two branches. The equality Sc{z] r) = implies that 
Im/ic(^;r) becomes zero. For Rez ^ this contradicts ljA.2|l . As before, the case 
Rez — must be checked separately. This can happen via (|A.3|1 . Altogether, we 
conclude that £c{z] r) is analytic and nonzero on each of its two branches. 

Before we can make statements about the eigenvalue A(fc,z), we have to pay 
attention to the eigenvector u(fc,z;r), which is the solution of (|21|l and H22|) . As 
is manifest from H17|l . the dependence of the operator L{s\y) on the real variable 
s is smooth, except for the points where /i(s; r) equals zero. The corresponding 
divergencies do not appear in u(fc, s; r), owing to the normalisation (|19() . This incites 
us to suppose that u(fc, s; r) is regular and smooth for all real s. Consequently, 
u(fc,z;r) is analytic in a certain strip around the real axis. Outside this strip the 
operator i(z;r) is nowhere singular, so it seems reasonable to assume that u(fc,z;r) 




(A.l) 



From this result the inequality 

|Imft,(z;r)| > |lm(z^) 



(A.2) 



h{ib\Y) < -wo(i-) 



(A.3) 
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is an entire function in the complex z plane. Then the same goes for [u(fc, z*; r)]*, by 
the rules of complex conjugation. 

The results for /ic(z;r) and £c{z\y)^ which were proved above, illustrate that the 
continuum limit does not create any new singularities outside the real axis. It merely 
modifies the character of the singularities that already exist on the real axis. Poles 
unite and form a branch cut. We assume that also for u(fc, z;r) the continuum limit 
does not create any new singularities. Said differently, the eigenvector remains an 
entire function. Note that all of our assertions on u(fc, z; r) are trivial for the case of 
a homogeneous dielectric, because then the dependence of u(fc, z; r) on z is absent. 

In H25I) the eigenvalue X{k, z) is expressed in terms of a spatial integral over the 
volume V . The above results tell us that the integrand is analytic in the whole complex 
z plane, except for the real axis. Since V is finite, the integral is convergent if the 
integrand is bounded. We therefore may conclude that \{k, z) is analytic as long as 
Imz differs from zero. The foregoing argument can be repeated for Ac(A:,z), so this 
function is analytic on both of its two branches. 

The identity [A(/c, z*)]* = \{k, z) allows us to put forward 

'\{k,z) - X{k,z*y 



ImA(A:, z) = Imz 



(A.4) 



Let us make Imz small, and choose Rez such that h{z:r) is nonzero for Imz = 0. We 
then avoid divergencies on the real axis. To elaborate l|A.4|l we invoke the identity 

^h^^ /"dr[u(fc,z*;r)]*-u(fc,z;r)= /" dr [u(fc, z*; r)]* • u(fc, z; r)^^^^^. (A.5) 

It can be proved with the help of 125|1 and partial integration. Before carrying out 
the replacement s — > z in H25|) . one has to multiply the left-hand side by the norm 
(u(fc, s), u(fc, s)). After use of (|A.5p the identity ljA.4|) attains the form 



ImA(fc,z) w Im(z2) 



dr |u(fc, z; r)| 



2 dz^e(z\ r) 



dz2 



(A.6) 



Im2=0 



where Imz is small and h(z\ r) must be nonzero for Imz = 0. 
The result (jA.6(l invites us to employ the inequality 

> 1 (A.7) 

which follows from (|16|l . This brings us to 

|ImA(fc,z)| > |lm(z2)| (A.8) 

where Imz is small and Rez is nonzero. Hence, the function A(fc, z) surely differs from 
zero in the vicinity of the real axis. Note that the case Rez = is unimportant, 
because and force A(/c, z) to be real and negative on the imaginary axis. By 
taking the continuum limit of A(fc, z), we do not create zeros for small Imz, because 
the lower bound (jA.8|) does not depend on N and A. Hence, close to the real axis 
Ac(A:,z) is nonzero. 

The above findings on analyticity and location of zeros permit us to make a 
preliminary statement. We may claim that outside the real axis A~^(fc,z) is analytic 
whenever A(/c, z) differs from zero. Therefore, in the vicinity of the real axis A~^(/c, z) 
is analytic. Now the question arises of what happens further out in the complex plane, 
and on the real axis itself. 
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For the case of a homogeneous dielectric X{k, z) has {2N + 4) zeros and {2N + 2) 
poles. AU of these are located on the real axis, as follows from (jA.4p . Under a smooth 
transition from a homogeneous to an inhomogeneous dielectric, the (4iV + 6) zeros and 
poles cannot leave the real axis. This has been demonstrated above. What we now 
assume is that the transition does not generate any new zeros or any new poles. Then 
X~^{k,z) is meromorphic, with all poles located on the real axis. Moreover, Ac(fc,z) 
will not possess any zeros on each of its two branches. 

To carry out a check on our assumption, we utilise the argument principle j24j . 
For A(fc, z) it can be formulated as 

^dM^^^^^ ^ (A.9) 

On the right-hand side the number of poles of A(fc, z) is subtracted from the number of 
zeros of A(fc, z). The contour C2 is a large circle centered around the origin. To prove 
()A.9(I we appeal to I) A. 5(1 once more. The ensuing integrand can be simplified with 
the help of IjB.ip and the fact that e(z; r) converges to unity for \z\ large. Then (|A.9|I 
boils down to the condition J^^dz/{2Triz) = 1. This elementary integral completes 
our consistency check. 

Our last job is the analytic continuation of hc{z). The symmetry (|67ll allows 
us to extend the contour of (|63|l to the negative real axis. By assumption, is 
regular and smooth for all real z, so it must be analytic inside a certain strip around 
the real axis. Keeping Imz positive as required by (|66|l . we shift in (|63|l the contour 
downwards, until it runs below the pole at lo = —z. Upon evaluating the residue we 
acquire an alternative form of l(63|l . given by 

If,, ^^/J^M 



The contour C4 is a straight line, running parallel to the real axis at Imuj = —7, with 
7 positive. 

Rather than the upper half plane, the strip —7 < Imz < 7 constitutes the region 
where the representation (|A.10|I is analytic. We thus have succeeded in finding an 
analytic continuation of ((63|l below the real axis. From ((A.10|l we infer the symmetry 
relation 

[K{-Z*t ^hrXz). (A.ll) 

Note that the operation z —t —z* maps the strip —7 < Imz < 7 onto itself, so in 
employing IjA.llll domain questions do not arise. In view of (|16|l and (|25|l . as well as 
the fact that u(fc, z) is an entire function, the analytic continuation of edij-i + iri) and 
Xc{k,Lo + iry) can be performed on the basis of (|A.10|I . 

All properties of h and related functions, which are needed in the main text, 
have now been proved or made plausible. In essence only one basic assumption on 
the partial differential equation ((22|l underlies our discussion. For the eigenvector 
u(fc,z;r), the eigenvalue A(fc,z), as well as its inverse A~^(fc,z), both the continuum 
limit, and the transition from a homogeneous to an inhomogeneous dielectric, do not 
give birth to singular points out in the complex plane. As a justification, we point out 
that all singular points of the operator L(z; r) lie on the real axis, regardless of the 
decision to carry out one of the two afore-mentioned procedures. Of course, the check 
Tg)) provides further support for our views. 
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Appendix B. Verification of two canonical commutation relations 



In subsection 3.2 we commenced the verification of the upper two commutators H41|l . 
Upon making the transition to complex integration, we arrived at the sufficient 
conditions (|15|l - l|T7|) . The reflection principle [\{k,z*)]* = X{k,z) ensures that the 
integrals (|45|) and 1)47(1 are invariant under complex conjugation, combined with the 
interchange of position vectors and tensorial indices. This invariance is imposed by 
the tensorial delta functions figuring in the right-hand sides of (|45ll and H47|) . We first 
get over with 145|) . which is the simplest condition. 



(a) 


Ci 








(b)_- 









(c) 






^ 












Figure Bl. Integration contours in the complex plane z. Horizontal axis: 
Im2 = 0; vertical axis: Rez = 0. (a) Contour C\ consists of {2N + 4) circles of 
infinitesimal radius; each circle runs in counterclockwise direction and is centered 
around a point of the set {z = il{k, I)}, (b) Contour C2 consists of a single circle 
of large radius, centered around the origin, (c) Contour C3 consists of two straight 
lines running above and below the real axis at infinitesimal distance. 

With the help of ((21|) and the transformation z — z, we can include the poles at 
z = — n(fc, I) in the contour. In front of the integral a factor of 1/2 appears. The new 
contour C2 (see figure) is a circle z = i?exp(z^),0 < 4> < 2it. All poles z = ±.^{k,l) 
lie in the interior of C2, because R can be taken arbitrarily large. Note that, except 
for the poles z = ±fl{k,l), the integrand is everywhere analytic. A crucial advantage 
of C2 over Ci is that one can benefit from asymptotic results. If \z\ is large, (|22|l has 
the algebraic form 

u{k, z;r) ^ u{k, z;r). (B.l) 



A(fc, z 

The condition (|45l) becomes 



E 

k 



f ^i[uT(fc,z*;r')]*UT(fc,z;r) = JT(r'-r). (B.2) 

Use of (|B.1)I has modified the analytical structure of the integrand. The poles at 
z = ±f2(fc, I) have given way to a single pole at the origin. 

We effectuate a last deformation of our integration contour. Let C3 enclose the 
real axis by means of two straight lines running from +00 + irj to —00 + irj, and from 
—00 — ir] to +00 — ir], where t] is infinitesimally positive. The integrand in l|B.2p is 
still analytic outside the real axis, so integration along C3 instead of C2 is permitted. 
Now we are in a position to appeal to the completeness relation (|20|l . This leaves us 
with the condition J dz/{2Triz) = 1, which is indeed true for the contour C3. 

To get to grips with the analytical structure of the integrand of , we make 
the replacement s ^ z in l|22(l . and write the result as 
Q!^(r)z^u(fc, z; r) 
eop(r)/i(z;r)A(fc, z) 
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\-^{k,z) {-c^V X [V X u(/c,z;r)] + z'^\i{k, z]r)] - u(fc,z;r). (B.3) 

Obviously, the function on the left-hand side has the same analytic structure as 
A~^(fc,z). Hence, (|46|) does not really cause any new complications. We first include 
the regular point z = in the contour Ci, and next perform the substitution ljB.3|l . 
Deformation of contours eventually leads to the condition / dz/[zh{z\Y)\ = 0, where 
Cs is the contour. As shown in appendix A, all poles of the meromorphic function 
/i~^(z;r) are located on the real axis. After deformation of C3 to the large circle C2, 
the integral indeed disappears, because h{z;r) behaves as z^ for j^l large. 

Since the integrand of (|47|l has two factors of h in the denominator, its analytic 
structure must be carefully investigated before any contour deformations can be 
undertaken. This can be done by adding a new contribution to the integrand, the 
analytic structure of which is determined by the function h^^{z;r'). Instead of the 
integrand itself, we consider the sum 

, . _ z'^[uik,z*;r')]*u{k,z;r) zhop{r)[u{k, z*;r')]*u{k, z;r) 
^ ' ' ' ^ h{z;r')hiz;r)\{k,z) ^ a^r)h{z;r') ' 

The reason is that the poles of the function f can be easily located. Upon employing 
ljB.3|l twice, one recognises that f has the same analytic structure as \~^{k,z). Thus 
the integrand of H47() has poles both for A(/c, z) — Q and for h{z;Y') — 0. Again, we 
include the regular point z = in the contour Ci. Now we can start shifting contours. 
At each point z = Q.{k, I) the second contribution to f is regular, so integration along 
Ci yields zero. Consequently, in H47|) the integrand may be exchanged for f/z. Next, 
one can replace Ci by C2, if the transformation z — ^ — z is used. Upon employing 
IIB.411 once more, one arrives at the condition 



k 



dz 



■Ic2 27i"i 



/i(z; r)A(fc, z) a^(r) 



z[u(fc,z*;rO]*u(fc,z;r) _ ^op(j).r, . .-r 



We apply (|B.1|) to the first term between square brackets. Subsequently, we 
replace C2 by C3, so that (|20|l comes into play again. Eventually, the integrals 
/ dzz/[h{z;r)h{z;r')] — and / (izz/[27ri/i(z; r)] = 1 must be proved. As before, 
this can be done by passing over from the contour C3 to the large circle C2. 
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